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Abstract. The evolution of young compact star clusters 
is studied using Y-body simulations in which both stellar 
evolution and physical collisions between stars are taken 
into account. The initial conditions are chosen to represent 
R136, a compact star cluster in the 30 Doradus region 
of the Large Magellanic Cloud. The present runs do not 
include the effects of primordial binaries. 

We find that physical collisions between stars in these 
models are frequent, and that the evolution of the most 
massive stars and the dynamical evolution of the cluster 
are closely coupled. In all cases, a single star grows steadily 
in mass through mergers with other stars, forming a very 
massive ( }t 1OOM ) star in less than 3-4 Myr. The growth 
rate of this runaway merger is much larger than estimates 
based on simple cross-section arguments, mainly because 
the star is typically found in the core and tends to form 
binaries with other massive stars there. The runaway is 
"rejuvenated" by each new collision, and its lifetime is 
extended considerably as a consequence. Observationally, 
such a star will appear in the Hertzsprung- Russell diagram 
as a blue straggler. When the runaway forms a black hole, 
the binary in which it is found is usually dissociated. 

We further investigate the sensitivity of the runaway 
to different formulations of mass loss from high-mass main 
sequence stars. We find that, while the runaway process 
is less pronounced in the presence of strong stellar winds, 
the basic effect persists even in the face of large mass loss. 

Key words: binaries: close — blue stragglers — stars: 
evolution — stars: mass losss — globular clusters: general 
— globular clusters: 30 Doradus 



Send offprint requests to: Simon Portegies Zwart: 
spz@grape.cu-tokyo.ac.jp 

* Japan Society for the Promotion of Science Fellow 



1. Introduction 

Physical collisions between stars are not rare in the central 
regions of star clusters or galaxies. In some star clusters, 
and in the cores of many galaxies, stellar collisions are 
likely to play an important role in the formation of exotic 
objects such as blue stragglers (Sanders 1970; McNamara 
& Sanders 1976), X-ray binaries (Fabian et al. 1975) and 
millisecond pulsars (Lyne et al. 1987; 1988). Collisions 
may also be responsible for the color gradient observed 
in several post-collapse globular clusters (Djorgovski et 
al. 1991), and the formation of central black holes in galac- 
tic nuclei (Quinlan & Shapiro 1990; Quinlan et al. 1995, 
Lee 1995). Collisions in proto-clusters may be responsible 
for the formation of massive stars (Bonnell et al. 1998), 
and possibly even the entire mass spectrum (Silk & Taka- 
hashi 1979; Allen & Bastien 1995; Price & Podsiadlovski 
1995) 

In order to quantify the effect of collisions on the 
evolution of stars and the corresponding changes in the 
stellar population, Portegies Zwart et al. (1997, here- 
after paper I) performed population synthesis calculations 
which included stellar collisions. In those calculations the 
stellar number density was held constant, thus excluding 
the possibility of any interplay between the dynamical evo- 
lution of the cluster and collisions between stars. 

In reality, both the stars and the parent cluster evolve 
on comparable time scales, and cluster dynamics and stel- 
lar evolution are quite closely coupled. For example, mas- 
sive stars tend to segregate to the core due to dynamical 
friction, increasing their collision probability. Their colli- 
sion products are even more massive, leading to the pos- 
sibility of runaway merging (Lee 1987; Quinlan & Shapiro 
1990), if the collision rates can remain high enough in 
the few Myr before the stars explode as supernovae. How- 
ever, these rates are determined by the dynamical state of 
the cluster, which is in turn strongly influenced by stel- 
lar mass loss. The only way to treat this intimate cou- 
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pling between stellar collisions and cluster dynamics is 
to perform JV-body simulations in which the stars are al- 
lowed to evolve and collide with one another in a fully 
self-consistent way. 

In this paper we report the results of a series of iV-body 
simulations modeling young and compact star clusters, 
such as R 136 in the 30 Doradus region in the Large Mag- 
ellanic Cloud. This cluster is particularly interesting be- 
cause a strong coupling between stellar evolution and stel- 
lar dynamics may exist. In addition to this, excellent ob- 
servational data is available. Many unusually bright and 
massive stars (e.g. Massey & Hunter 1998) are present in 
R136 which, due to the high central density of 10 6 to 10 7 
stars pc -3 , are likely to interact strongly with each other. 

The numerical method is discussed in Sect. 2. Sect. 3 
describes in more detail the initial conditions for our mod- 
els. In Sect. 4 the results are presented; they are discussed 
in Sect. 5. Briefly, we find that runaway collisions of mas- 
sive stars can occur. The most massive star grows in mass 
through merging with other stars until it collapses to a 
black hole. The growth rate of this star is much larger 
than estimates based on simple cross-section arguments, 
because the star is typically found in the cluster core, and 
tends to form binaries with other massive stars. 

2. Numerical method 

The A^-body integration algorithm, used in this paper, is 
described in Sect. 2.1. In Sect. 2.2 we describe how the 
evolution of stars are calculated; the effect of collisions on 
the evolution of stars is described in Sect. 2.3. 

2.1. The N-body integrator 

The iV-body portion of the simulations is carried out using 
the kira integrator, operating within the Starlab software 
environment (McMillan & Hut 1996; Portegies Zwart et al. 
1998). Time integration of stellar orbits is accomplished 
using a fourth-order Hermite scheme (Makino & Aarscth 
1992). Kira also incorporates block timesteps (McMillan 
1986a; 1986b; Makino 1991) special treatment of close 
two-body and multiple encounters of arbitrary complex- 
ity, and a robust treatment of stellar and binary evolu- 
tion and stellar collisions (see below). The special-purpose 
GRAPE-4 (Makino et al. 1997) system is used to acceler- 
ate the computation of gravitational forces between stars. 
The treatment of stellar mass loss is as described in Porte- 
gies Zwart et al. (1998). A more complete description the 
Starlab environment is in preparation. 

2.2. Stellar evolution 

The evolution of stars is taken from the prescription by 
Portegies Zwart & Verbunt (1996, Sect. 2.1). However, 
some changes are made to the mass loss in the main- 
sequence stage for massive stars. 



2.2.1. Mass loss from main-sequence stars 

The original equations from Egglcton et al. (1989), on 
which the stellar evolution model is based, ignore mass 
loss during the main sequence stage. However, for stars 
more massive than 25 M© , mass loss on the main sequence 
can be substantial. We use three different prescriptions to 
investigate the effect of main-sequence mass loss. 

The first prescription simply follows Egglcton ct al., 
and no mass is lost on the main-sequence. In these cases, 
a massive star loses its entire hydrogen envelope when it 
leaves the main-sequence and becomes a Wolf-Rayet star. 
We refer to this prescription as no mass loss. 

In the second prescription the mass loss rate for a mas- 
sive main-sequence star is taken to be constant in time, in 
such a way that, as it leaves the main sequence, the star 
has lost its entire hydrogen envelope. We refer to this type 
of mass loss as constant mass loss. 

In the third and most realistic treatment, a massive 
star loses its hydrogen envelope during the main-sequence 
phase according to the law 

moc-t 1 ' 5 . (1) 

This treatment is supported by model computations for 
massive stars by Schaerer et al. (1999). We refer to this 
type of mass loss as moderate mass loss. 

In several cases the mass of collision products exceeds 
120 M© (the most massive evolutionary tracks available, 
see Schaerer et al. 1992). Very little is known about the 
evolution of such massive stars (see Stothers et al. 1997; de 
Koter et al. 1998 and Figer et al. 1998). We assume that for 
such high mass stars the lifetime and the radius depends 
weakly on mass (see Langer et al. 1994). In our model a 
100 M© star has a main-sequence lifetime of 3.08 Myr; a 
150 M Q star lives for 2.98 Myr. 

2.2.2. Supernovae and velocity kicks 

A star with a mass larger than 40 M Q leaves a black hole 
after ejecting its envelope during the main-sequence and 
Wolf-Rayet phases. The mass of the black hole is com- 
puted as mbh = 0.35m — 12 M©, where mo is the initial 
mass of the star. For a star whose mass increases due to 
collisions, to is the highest mass reached by the star. 

Stars with masses between 8M© and 40 M Q become 
neutron stars. At birth a neutron star receives a high 
velocity 'kick' in a random direction. The magnitude of 
the velocity kick is chosen randomly from the distribution 
proposed by Hartman (1997). This distribution is flat at 
velocities below 250 km s~ , but has a tail extending to 
several thousand km s _1 . 

Stars with masses less than 8 M© become white dwarfs. 
The mass of the white dwarf equals the core mass of its 
progenitor at the tip of the asymptotic giant branch. 
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2.3. Stellar collisions 

A collision is assumed to occur when two stars (i and j) 
approach each other within a distance d = 2(r,+r 7 ), where 
ri and Tj are the radii of the stars involved. 

2.3.1. The collision product 

A detailed description of the treatment for collisions is 
given in paper I (section 3.3). Here we summarize the pre- 
scription for collisions between main-sequence stars. 

A collision between two main-sequence stars with 
masses nij and nij results in a single rejuvenated main- 
sequence star with mass m, + rrij . Smooth-particle hydro- 
dynamic simulations of collisions between main-sequence 
stars indicate that at maximum a few percent of the total 
mass is lost (see e.g.: Lai et al. 1993; Lombardi et al. 1995; 
1996). Consequently, mass loss during the merger event is 
ignored. 

The collision results in a reduction of the age of the col- 
lision product. The age reduction factor / rod is computed 
from the mass of the most massive of the two colliding 
stars rrii and the mass of the collision product (Meurs & 
van den Heuvel 1989) 



/ re d(TOi,m.j) 



mt + rrij T ms (mi) 



(2) 



Here T ms (mi) is the main-sequence lifetime of a star with 
mass rrii. The new age of the collision product is computed 
with i*(m; + rrij) = / re d(^i, rnj)U(mi). 

As an example, suppose that, at t = 5 Myr, a star 
with rrii = 20 M collides with a second star with rrij = 
8M Q . Both stars lie on the main sequence, and both are 
experiencing a collision for the first time (i.e. i+ = 5 Myr 
for both stars). For a 20 M star, T ms (2OM ) ~ 7.7Myr. 
The collision results in a main-sequence star with a mass 
of 28 M , for which r ms (28M ) ~ 5.4Myr. The new age 
of the collision product is computed using Eq. (2) , which 
in this case gives i*(28M ) ~ 2.5 Myr. 

3. Selection of initial conditions 

We selected the initial parameters for the models to mimic 
class of star clusters similar to the Galactic cluster NGC 
3606 or the young globular cluster NGC 2070 (R136) in 
the 30 Doradus region of the Large Magellanic Cloud. 

3.1. The star cluster R 136 in the 30 Doradus region 

The half-mass radius (rh m ) of R136 is about 1 parsec 
(Brandl et al. 1996), and the core radius is <~ 0.02 pc 
(Hunter et al. 1995). The total mass M ~ 2 x 10 4 M . 
With an assumed mean mass of 0.6 M the cluster thus 
contains about 35 000 stars. The corresponding central 
density is of the order of lO 6 M pc- 3 . The age of R 136 
is - 3-4 Myr (Campbell et al. 1992). The Galactic star 



cluster NGC 3606 is somewhat smaller in size and its to- 
tal mass is larger resulting in a denser core (Moffat et al. 
1994; Drissen et al. 1995). 

For both clusters, the tidal effect of their parent galaxy 
(for R 136 that is the Large Magellanic Cloud) is small. 
Assuming that the mass of the LMC is Mlmc ~ 10 9 M 
and the distance from the center of the LMC is ~ 1 kpc, 
the tidal radius, 



n 



( M \ 1/3 



(3) 



is ~ 20 pc, much larger than the half-mass radius of either 
cluster, justifying our neglect of tidal effects. 

Calculations are performed using 12k and 6k stars. 
Therefore we need to scale the dynamical timescale and 
the collision cross section to mimic the evolution of a star 
cluster with larger N. This scaling is discussed in the fol- 
lowing two sections. 

3.2. Scaling the dynamical timescale 

The evolution of an isolated star cluster is driven by two- 
body relaxation. Therefore, we set up the initial model 
so that it has the same relaxation timescale as the real 
cluster. The relaxation time is calculated with 



i rlx cx 



N 



logiofrW) 



thu 



(4) 



Here 7 is a scaling factor, introduced to model the effects 
of the cut-off in the long range Coulomb logarithm (see 
Giertz & Heggie 1996; 1994). Here ihm is the half-mass 
crossing time of the cluster is 



(5) 



Here rhm is its half mass radius. 

The radius of the scaled model r moc i i is then computed 
by substitution of Eq. (4) into Eq. (5): 



^*modcl 



Wrcal A' 7 ' / ln( 7 A rcal ) V 2/3 



model 



model , 



^real* 



(6) 



Here r roa i is the radius of the real cluster. We decided to 
use 7 = 1. Note that Eq. (6) is relatively insensitive to the 
exact choise of 7. 

3.3. Scaling the collision cross section 

The model clusters should have the same collision rate 
per star as the real cluster. Scaling the initial conditions to 
assure that the model cluster has the same relaxation time 
causes it to be larger than the real system (Eq. 6) . The 
correct collision rate per star is then obtained by scaling 
the sizes of the stars themselves. 
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by 



"coil oc n c av. 



(7) 



Here n c is the number density of the stars in the core, 
a is the collision cross section (for approach within some 
distance d), and v is the velocity dispersion. These are 
given by the following proportionalities: 



N 

n c oc — , 

ri. 



a oc d 2 H — 7T. 



(8) 



Here r c is the cluster's core radius. We will neglect the 
d 2 term in the cross section. Expressed in real units and 
assuming scaling according to Eq. (6), we may write 



N 



oc N 2 , 
v oc TV 2 / 3 . 



The number of collisions then becomes 
ncoii = dN A '\ 



(9) 



(10) 



The distance at which a collision occurs therefore 
scales as 



d oc N- A '\ 



(11) 



The names used to identify our models are defined as 
follows: We start with the number of stars: 12k for mod- 
els with 12288 stars, and 6k for models with 6144 stars. 
The next integer identifies the selected value for Wo- The 
next letter indicates the stellar mass loss model: 'A' for no 
mass loss, 'B' for constant mass loss and 'C for moder- 
ate mass loss (see Sect. 2.2.1). The final number gives the 
initial half-mass relaxation time in millions of years. Two 
12k6C10 models were computed. To distinguish between 
them, the second is identified as 12k6C10'. All computa- 
tions were continued until t = lOMyr. 
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Fig. 1. Initial parameters for the 6k runs. 



3.4- The models 

We performed 4 runs with 12k stars and 7 runs with 6k 
stars. All 12k runs start from the same initial conditions, 
except for the treatment of mass loss on the main se- 
quence. For 6k models we also change the initial relaxation 
time and the initial density distribution. Table 1 and Fig. 1 
summarize the initial conditions. 

All simulations start at t — by assigning masses of 
stars between 0.1 M© and 100 M© from the mass function 
suggested for the Solar neighborhood by Scalo (1986). At 
the high-mass end this mass function is rather steep; 



N{m) oc m 



-2.82 



(12) 



and the mass function turns over at around 0.3 M©. The 
median mass of this mass function is about 0.3 M©, and 
the mean mass is about 0.6 M©. We generate the mass dis- 
tribution using the random sampling technique suggested 
by de la Fuente Marcos et al. 1997). 

The initial density profile and velocity dispersion for 
the models with 12k stars are taken from a King (1966) 
model with Wo = 6. We chose r v i r = 0.31 pc, which results 
in a core radius r COIC ~ 0.072 pc and a core density of 
/Wc ~ 3.6 x 10 5 M© pc -3 . The central velocity dispersion 
for these models is about 8.7 km s _1 and the initial half- 
mass relaxation time t T i K ~ lOMyr. 



4. Results 

In this section we describe our results. Stellar collisions 
and the growth of massive stars are discussed in Sect. 4.1; 
Sect. 4.2 discusses evolution of the cluster structure, while 
Sect. 4.3 gives a detailed description of the evolution of 
the runaway collision product in one particular model 
(12k6A10). The results of the models with 6k stars are 
described in Sect. 4.4. 



Table 2. Results on the N = 12 k runs. The first column gives 
the model name. The following columns give the time of the 
first collision, the total number of collisions in the computation, 
the number of collisions in which the most massive star partic- 
ipates (A r M max ) and the number of supernovae which occurred 
during the computations. 



Model 


^l at coll 

[Myr] 


iVcoll 


^M max 


N Bn 


12k6A10 


1.7 


11 


11 


5 


12k6B10 


0.9 


10 


5 


2 


12k6C10 


0.3 


21 


15 


4 


12k6C10' 


1.1 


15 


10 


9 
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Table 1. Overview of initial conditions for the runs with N = 12k (first four rows) and N = 6k (last 7 rows). The first 
column gives the model name, starting with the number of stars: '12k' for the 12k models and '6k' for the 6k models. The 
following columns give the number of stars, the initial dimensionless King Wo, the relaxation time t T i x and the crossing time at 
the half-mass radius ihm (both in millions of years). The next columns give the initial core radius r CO rc, the half-mass radius 
fhm (both in parsec) and the average core density (in Mo/pc 3 ). The last three columns give the number of stars initially more 
massive than 19 M Q (the turn-off mass for a 10 Myr old star cluster) and 25 M (the mass limit for a Wolf-Rayet star) , and the 
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0.07 


0.25 


5.56 


7 


4 
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6 
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6144 
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{.1. Runaway merging 
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Fig. 2. Diagram of the encounter history for the most massive 
stars in models 12k6A10 (dotted line), 12k6B10 (dash dotted 
line), and models 12k6C10 (solid line) and 12k6C10' (dashed 
line). The horizontal axis gives time in millions of years, the 
vertical axis shows the mass of the collision product. The "•" 
symbol indicates the moment a star leaves the main sequence 
and becomes a Wolf-Rayet star. The "*" indicates the moment 
the star collapses to a black hole. The vertical line after the 
* shows the mass lost in the supernova. The evolution of the 
star after the supernova is not shown. Each collision is indi- 
cated with a small o. For some models two stars are initially 
displayed. In these cases, the two stars ultimately merge, but 
only after each experiences collisions with less massive stars. 
For example, in model 12k6C10, the 62 Mq star experiences 
a collision with a lower mass star before it encounters the 
70 M Q star, which itself had experienced 4 collisions prior to 
the merger. 



Figure 2 gives the evolution of the mass of the most mas- 
sive star for all TV = 12 k runs. In all cases, more than 
10 collisions occurred during the first 5 Myr. In model 
12k6A10 (no mass loss), a total of 11 collisions resulted in 
a star with a mass of 182 M Q at the point when it left the 
main sequence. 

From Fig. 2 we can see that the lifetime of the most 
massive star is considerably larger than its "natural" 
main-sequence lifetime of about 3 Myr. This extension of 
the lifetime is caused by rejuvenation through merging. 

In other models, two massive stars collide with rel- 
atively low-mass stars before they eventually find each 
other. In models 12k6C10 and 12k6C10' this happens at 
t ~ 1.7 Myr, in model 12k6B10 at t ~ 3 Myr. 

Table 2 gives information about collisions in the N — 
12 k models. It typically requires about 1 million years - 
the time needed for the most massive stars to segregate 
to the core- before the first collision occurs (see Table 2). 
These most massive stars participate in more than 70% 
of all collisions. The most massive star in model 12k6C10 
experiences its first collision at an earlier epoch because 
that star happened to be born in the core. 

The number of collisions occurring in these runs is 
far larger than simple theoretical predictions. The rate 
at which stars in a cluster experiences collisions can be 
estimated as n 2 c {av) (Spitzer 1987). Here a is the collision 
cross section, v is the velocity dispersion in the cluster and 
n c is the number density of stars. Following the derivation 
in paper I by adopting a Maxwellian velocity distribution 
with velocity dispersion (v) and cross section a cx md/v 2 , 
the number of collisions is the cluster per Myr is expressed 
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as (see paper I, Eq. 14) 



0.36 



(. 



V10 4 M Q pc- 3 
d \ /km s~ 



pc 



2 (to) 



[Myr" 



(13) 



Here (to) is the mean stellar mass. 

For the N = 12 k runs this results in about 0.3 colli- 
sion/Myr, or about 3 collisions during the entire simula- 
tion, assuming that the cluster parameters do not change 
in time. As will be discussed in Sect. 4.2 below, the core 
density in fact drops by about an order of magnitude dur- 
ing the first 4 million years. If we take this effect into ac- 
count, the expected number of collisions is less than unity. 
The actual number of collisions in the 12k simulations ex- 
ceeds 10. The major cause of this large discrepancy is mass 
segregation, which concentrates massive stars in the core. 

The importance of mass segregation is illustrated in 
Figs. 3 and 4. Figure 3 gives the theoretical probability dis- 
tribution for stars with mass TO, pr i m to collide with lower- 
mass stars of mass to scc , for the initial Scalo (1986) stel- 
lar mass distribution. Figure 4 presents the distribution of 
collisions actually observed in our simulations. 




0. 1 



10 



M . [Mr,] 
przin L J 



Fig. 3. Expected relative collision rate between a primary and 
a secondary star for the adopted mass function and mass-radius 
relation for main-sequence stars. The velocity dispersion is 
taken the same for all masses and Eq. (13) is used to compute 
the relative encounter probabilities. The shading in linear in 
the encounter probability, with darker shades for higher prob- 
abilities. Since the probability distribution is symmetric line 
the axis of equal mass, only the lower half is displayed. 



The differences between Fig. 3 and Fig. 4 are striking. 
In the A-body simulations, massive stars completely dom- 
inate the collision rate, while theory predicts that the ma- 
jority of collisions should occur between stars of relatively 
low mass (~ 0.7 M Q ). 

Figure 5 shows the effect of merging on the statistics 
of supernovae. The total number of supernovae is reduced 



100 r- 



1 r- 




0. 1 



Fig. 4. Primary and secondary masses of colliding stars for 
all simulations. A total of 125 collisions are combined in this 
graph. Note that the axes extend to 200 M© , compared with 
10 M in Fig 3. 



due to merging of massive stars, and the explosions are 
delayed relative to expectations because of rejuvenation. 



T [Myr] 



Fig. 5. Upper symbols (•): Supernova history for model 
12k6A10 over the first 16 Myr of the evolution of the stellar 
system. Lower symbols (o): the expected (scheduled) super- 
novae if dynamics did not play a role. 



4-2. Evolution of cluster structure 

The core and the cluster as a whole expand with time 
(see Fig. 6; for technical reasons, fewer snapshots of the 
first model were stored, leading to lower temporal resolu- 
tion in the data displayed here). For model 12k6A10 this 
expansion is almost completely driven by binary heating. 
For the other models, mass loss in the stellar winds of the 
massive stars also drives the expansion. Although in these 
models the cluster loses a modest 4% of mass in the first 
10 Myr, this mass is lost from deep inside the potential 
well of the cluster and affects the dynamics significantly. 

Figure 7 shows the evolution of the mean stellar mass 
in the core. The initial mean mass in the cluster is (to) « 
0.6 M q . In the core, (to) increases to about 1.1 M Q in 
1 Myr due to mass segregation. In model 12k6AlO, (m) coro 
keeps increasing until a maximum is reached at t ~ 5 Myr. 
In model 12k6C10, the maximum is reached at t ~ 2 Myr. 
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Fig. 6. Core and Lagrangian radii (in parsec) as functions of 
time (in Myr) for models 12k6A10 (dotted lines) and 12k6C10 
(solid and dashed lines). The lowest lines represent the core 
radius, the higher lines give the 25%, 50% and 75% Lagrangian 
radii. 

This difference is in part due to the different treatments 
of mass loss from main-sequence stars. 
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Fig. 7. Mean stellar mass in the core as a function of time for 
models 12k6A10 (dots) and 12k6C10 (solid). The mean mass 
in the cluster is about 0.61 M , and decreases by only a few 
percent over the course of the simulation. 

Figure 8 shows the evolution of the core density p COIC 
for models 12k6A10 and 12k6C10. The expansion of the 
core causes /9 CO rc to decrease during the first 4 million 
years. 

4-3. Anatomy of a collision sequence 

The first binaries in this model are formed shortly after 
the start of the simulation. One of these binaries is formed 
from the most massive star, of mass 57 M Q , and a compan- 
ion of 26 M , with an initial semi-major axis of 0.03 pc. 
After more than ten collisions, the 57 M© star grows to 
182 M Q , before turning into a black hole of 133 M©. 



Fig. 8. Central density (in M pc 3 ) as a function of time for 
models 12k6A10 (dotted line) and 12k6C10 (solid line). 

As an illustration of the merger process, Fig. 9 depicts 
a schematic reaction network of the runaway merger in 
model 12k6AlO. 

At about 1 Myr the binary containing the most mas- 
sive star encounters another binary. This second binary is 
dissociated, and its component stars are ejected from the 
cluster. The effect of this is noticeably in Fig. 10 in a steep 
rise in the binding energy of the remaining binary. Sub- 
sequently, a series of encounters with single stars results 
in three collisions involving the most massive star. The 
runaway collision product, still the member of a binary, 
devours two other stars between t = 2 Myr and 3 Myr. 

In the meantime a 27 M© star has formed a binary with 
another star. In the core an encounter between the two bi- 
naries results in a collision between the runaway merger 
and its 13 M© companion. This encounter dissipates most 
of the binding energy in the binary (see Fig. 10). The 
27 M© star is ejected from the core to return again as 
the member of a binary. 

The 27 M© star encounters the collision product again, 
the new companion of the 27 M© star collides with the run- 
away merger, and in addition the runaway merger collides 
again with its own companion; the 27 M© star takes its 
place. 

Slightly after t ~ 5 Myr the runaway merger collapses 
to a black hole after having consumed another 1 M© star. 
The binary survives but shortly afterward its companion 
(the 27 M© star) explodes in a supernova. The neutron 
star (remnant of the 27 M© star) receives a high velocity 
kick of ~ 100 km s _1 , dissociating the binary. The black 
hole (remnant of the runaway merger) is also ejected. 

Figure 10 plots the binding energy of the binaries and 
indicates some other important events for model 12k6A10. 
Similar figures are presented for model 12k6C10 in Fig. 11 
and for model 12k6A10' in Fig. 12 

Once binaries form, they gradually become harder un- 
til a maximum binding energy is reached. At that point 
the primary coalesces with its companion, removing its 
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Fig. 9. Reaction network from model 12k6A10, showing the evolution of the most tightly bound binary, which also contains 
the massive runaway merger. Time runs from left to right. Each line represents a star; the numbers attached to the lines give 
the masses of the stars in Mq. When two lines join (marked by circles), a binary is formed. A filled circle indicates that a 
strong interaction with another binary (not further specified) takes place. The masses of the binary components, rounded to 
Mq , are indicated within braces. The difference between the number of lines entering and leaving a circle indicates the number 
of collisions occurring. The 27 Mq star returns repeatedly to encounter the binary which contains the runaway merger. The 
"bh" slightly after t = 5 Myr indicates the collapse of the runaway merger to a black hole. The binary is finally dissociated by 
a supernova which ejects a neutron star ns (the former 27 M star) and a black hole (the runaway merger) . 
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Fig. 10. Energy stored in hard (Et > 5fcT, with kT the mean 
stellar kinetic energy) binaries as a function of time for model 
12k6A10. The numbers indicate the number of hard binaries at 
that instant; they are printed only at the moment they change. 
Small open circles at the top of the figure indicate the moments 
at which collisions occur. The larger symbols above them in- 
dicate the start of a Wolf-Rayet phase "•" , the formation of a 
black hole "*",and the formation of a neutron star (open star); 
a "o" indicates when a star starts its ascent of the super giant 
branch. 



Fig. 11. Energy stored in hard (Et > 5kT) binaries as a func- 
tion of time for model 12k6C10. See Fig. 10 for an explanation 
of the symbols. 



Most collisions occur between a member of a hard bi- 
nary and an incoming star. Following the collision, the bi- 
nary becomes softer. This is most clearly visible in Figs. 10 
and 12. During an episode without collisions, the binding 
energy of binaries rises at a rate of ~ 300fcT (about 0.2% 
of the binding energy of the cluster) per million years. 



binding energy from the system. The single star remain- 
ing after the merger captures a new companion, and the 
process repeats itself. Note that the term "binary" may 
be somewhat misleading here, as the runaway merger is 
usually the primary of a multiple system. The 8 "binaries" 
at t = 3.5 Myr in model 12k6A10 are in fact a hierarchical 
system of 7 lower mass companions orbiting the runaway, 
which at that instant has a mass of 151 M . 



4-4- Results of the 6k models 

Table 3 provides information on the calculations with 6k 
stars. The average number of collisions for the 6k runs 
with an initial relaxation time of 10 Myr is 7.0 ± 2.4. For 
the 12k runs the average number of collisions is 14 ± 4. 
The fact that the collision rate per star per unit time is 
about the same in both sets of models suggests that we 
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Fig. 12. Energy stored in hard (Et > 5kT) binaries as a func- 
tion of time for model 12k6C10'. See Fig. 10 for an explanation 
of the symbols. 



may be able to extrapolate our results to larger numbers 
of stars. 



Table 3. Overview of the 6k runs. The first five rows give 
information on the computations performed with no mass loss; 
the remaining two are from moderate mass loss runs. The first 
column gives the model name. Subsequent columns give the 
time of the first collision, the total number of collisions, and 
the number of collisions in which the runaway merging star 
participates (A^M max ). The final column gives the number of 
stars experiencing a supernova during the computation. 



Model 


^l Bt coll 

[Myr] 


Acoll 




AU 


6k6A5 


1.2 


21 


21 


1 


6k3A10 


2.6 


4 


1 


3 


6k6A10 


2.1 


9 


7 


1 


6k9A10 


0.6 


6 


4 


6 


6k6A20 


5.2 


3 








6k6C10 


1.2 


9 


5 


1 


6k6C20 


2.7 


5 


3 


1 



The results for models with different initial relax- 
ation times indicate that the collision rate is indeed in- 
versely proportional to the relaxation time, consistent 
with Eq. (13). On the other hand, the initial central den- 
sity has a rather small effect on the total number of colli- 
sions, even though the densities range over two orders of 
magnitude. 

Figure 13 shows a picture of the simulated cluster 
(model 12W6A10) from a distance of 2 parsec. An ani- 
mation of the cluster can be seen at the following address: 
http://www.sns.ias.edu/~starlab/research/30Doradu4Ke time needed for the first collisions to occur (about 
(http stands for uniform resource locator. The mpeg (Mov- 1 Myr) suggests that about 20 collisions could have oc- 
ing Picture Expert Group) animation shows the evolution curred in this cluster. The result of these collisions should 
of star cluster 12k6AlO from birth to about 7 Myr. In the therefore still be visible, possibly in the form of massive 



Fig. 13. Picture of 

one of the simulated clusters (model 12W6A10) from a dis- 
tance of 2 parsec. An animation of the cluster can be found 
at http://www.sns.ias.edu/~starlab/research/ or at the 
mirror http: //www . astro .uva.nl/~spz/30Doradus .html. 
The colors of the stars represent temperature. A collision prod- 
uct flares up in bright white. A supernova produces a bright 
violet star together with a slowly expanding shell-like struc- 
ture. 



beginning we look at the density center of the cluster from 
a distance of 30 pc (~ 100 r^c)- Then we zoom in with a 
velocity of about 29 km s _1 to a distance of 3 pc (~ 10rh c ). 

5. Discussion 

We have studied collisions in young star clusters with 6144 
stars and 12288 stars. The initial conditions of the models 
are selected to mimic stellar systems with a larger number 
of stars. The observed number of collisions is proportional 
to the number of stars, implying that our choice of scaling 
for time and stellar radius are appropriate. We therefore 
expect that the collision rate in our models can be extrap- 
olated to richer star clusters. The star cluster R136, for 
example, contains about 3 times as many stars as the mod- 
eled clusters. Scaling our results means that the collision 
rate in R136 is about 8 Myr" 1 . 

The absence of primordial binaries in our calculations 
possibly affects the collision rate and the collision counter- 
parts significantly (see Kroupa 1997; 1998). It is, however, 
not trivial to estimate the effect of the presence of a large 
fraction of primordial binaries; apparently it is not even 
easy to estimate the collision rate of a population of single 
stars correctly. 

5.1. The blue stragglers in R 136 

The small age of the star cluster R136 of ~ 3-4 Myr and 
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blue stragglers (Sanelage 1953; Leonard & Duncan 1990; 
Leonard 1995; Sills et al. 1997) in the cluster core. 

The three most massive stars in R136 have spectral 
type WN4.5 and appear to be younger than the other 
stars. Their age is estimated to be about 1 Myr (Massey & 
Hunter 1998; de Koter et al. 1997). These stars show violet 
absorption edges, which are common for late type (WN8 
and later) stars but highly unusual for these early types 
(Conti et al. 1983). Also striking is that these stars are 
unusually hydrogen rich (Massey & Hunter 1998). They 
are about an order of magnitude brighter than normal for 
such stars. Estimates for their masses range from 112 to 
155 M (Chlebowski & Garmany 1991; Vacca et al. 1996). 
Two of them lie well inside the core of the cluster; the third 
is at a projected distance of about 0.6 pc from the core. 

For a star cluster with an age of 4 Myr these three mas- 
sive stars appear as blue stragglers. It is therefore sugges- 
tive to identify the three most massive stars in the star 
cluster R136 as collision products. 



5.2. Black holes in dense star clusters 

When the runaway merger collapses to a black hole it is 
typically a member of a rather close binary. Upon dissoci- 
ation of the binary, the black hole is ejected from the core. 
Since the compact object is still considerably more mas- 
sive than average, mass segregation brings it back in the 
core within a few crossing times (sec e.g. Hut, McMillan, 
& Romani 1992). New close binaries can be formed once 
the black hole has returned to the core of the star cluster. 
After an episode of hardening the binary may become vis- 
ible as an X-ray source when the companion star starts to 
transfer mass to the black hole. Such a high-mass X-ray 
binary should be easily observable by X-ray satellites. The 
age at which such a binary can form is at least ~ 4 Myr, 
the minimum time needed for a black hole to form. It is 
likely to take considerably longer because the black hole 
has to return to the core after its ejection. 

The star cluster R136 is therefore "too young" for 
such a binary to exist. The star Mk34 at a distance of 
about 2.5 pc from the center of R136, however, is asso- 
ciated with a persistent X-ray source with a luminosity 
of ~ 10 36 ergs _1 (Wang 1995). Wang suggests that the 
binary contains a black hole of between 2.4 and 15 M Q 
accreting from the dense wind of its spectral type WN4.5 
Wolf-Rayet companion. This star, Mk34, can be classi- 
fied as a "blue straggler," as its estimated age is about 
lMyr, considerably smaller than the age of the cluster 
(De Marchi et al. 1993). 

Because R 136 is too young for such an X-ray binary 
to be formed from two collision products, it most likely 
formed from a primordial binary ejected from the cluster 
core following the supernova which formed the black hole. 



5.3. Collision rate 

The collision rates in our models are more than 10 times 
higher than simple estimates based on cross sections. In 
the computations with 12k stars, 14±4 collisions occurred 
in a timespan of about 4 Myr whereas only ~ 1 is expected. 

Furthermore, the cross section arguments imply that 
low mass (m £ 1M ) stars are most likely to collide. In 
our simulations, however, high-mass stars predominantly 
participated in encounters. The most massive star partic- 
ipates in numerous collisions with other stars. Typically, 
the mass of this runaway grows to exceed 120 M©. The 
rejuvenation of the runaway merger delays its collapse to 
a compact object following a supernova. Such a star could 
be visible in the core of young star clusters with a high 
density (IS <x blue straggler. 

The reason for the discrepancy between the formal 
cross-section arguments and the results of our simulations 
is the neglect of mass segregation and binary formation in 
the former estimates. In the simulations the most massive 
stars sink to the core due to dynamical friction within a 
few half-mass crossing times, and form close binaries by 
3-body interactions. The larger cross section of these bi- 
naries increases the collision rate and makes them favored 
candidates for encounters. 
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